function [CI_low, CI_high] = get_ci (est_main, est_ss, alpha, rate, N, b, method)

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Note: alpha needs to be specified as integers, e.g. [1; 5; 10]
    %
    % input:
    % est_main          ...main estimate theta
    % est_ss            ...sub-sample estimates
    % alpha             ...vector: percentiles of interest         
    % rate              ...convergence rate of the estimator 
    % N                 ...sample size
    % b                 ...block size
    % method            ...how should CI be computed?
    %
    % output:
    % CI_low            ...column vector: lower CI
    % CI_high           ...column vector: upper CI
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

factor = (b/N)^rate;
dim = size(alpha,1);
ones_dim = ones(dim,1);

switch lower(method)
   
   case {'backward','backwards'}

       perc_ss_low =  percentiles(est_ss, 0.5*alpha);
       perc_ss_high = percentiles(est_ss, 100*ones_dim - (0.5*alpha));
       
       CI_low = ones_dim*est_main + factor *  (perc_ss_low  - ones_dim*est_main);
       CI_high = ones_dim*est_main + factor * (perc_ss_high - ones_dim*est_main);


   case {'equal tailed','equal'}
     
       perc_ss_low =  percentiles(est_ss, 0.5*alpha);
       perc_ss_high = percentiles(est_ss, 100*ones_dim - (0.5*alpha));
       
       CI_low = ones_dim*est_main - factor * (perc_ss_high - ones_dim*est_main);
       CI_high = ones_dim*est_main - factor * (perc_ss_low  - ones_dim*est_main);
   
   case{'symmetrical','symmetric'}
     
       reps = size(est_ss,1);       
       T = abs(est_main*ones(reps,1) - est_ss);
       T_crit = percentiles(T, 100*ones_dim - alpha);
       
       CI_low = ones_dim * est_main - factor * T_crit;
       CI_high = ones_dim * est_main + factor * T_crit;
     
   otherwise
   
       error('Please specify CI method');

end
